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Introduction to Quantum Monte Carlo simulations for fermionic systems 
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We tutorially review the determinantal Quantum Monte Carlo method for fermionic systems, 
using the Hubbard model as a case study. Starting with the basic ingredients of Monte Carlo 
simulations for classical systems, we introduce aspects such as importance sampling, sources of 
errors, and finite-size scaling analyses. We then set up the preliminary steps to prepare for the 
simulations, showing that they are actually carried out by sampling discrete Hubbard-Stratonovich 
auxiliary fields. In this process the Green's function emerges as a fundamental tool, since it is 
used in the updating process, and, at the same time, it is directly related to the quantities probing 
magnetic, charge, metallic, and superconducting behaviours. We also discuss the as yet unresolved 
'minus-sign problem', and two ways to stabilize the algorithm at low temperatures. 

PACS numbers: 71.27.-(-a, 71.10.-w, 



I. INTRODUCTION 

In dealing with systems of many interacting fermions, 
one is generally interested in their collective properties, 
which are suitably described within the framework of Sta- 
tistical Mechanics. Unlike insulating magnets, in which 
the spin degrees of freedom are singled out, the interplay 
between charge and spin is responsible for a wealth of 
interesting phenomena (orbital degrees of freedom may 
also be included, but they add enormously to complexity, 
and shall not be considered here). Typical questions one 
asks about a system are related to its magnetic state (Is it 
magnetic? If so, what is the arrangement?), to its charge 
distribution, and to whether it is insulating, metallic, or 
superconductor. 

A deeper understanding of the interplay between spin 
and charge degrees of freedom can be achieved through 
models, which, while capturing the basic physical mecha- 
nisms responsible for the observed behaviour, should be 
simple enough to allow calculations of quantities compa- 
rable with experiments. The simplest model describing 
interacting fermions on a lattice is the single band Hub- 
bard modelQ, defined by the grand-canonical Haniilto- 



where t is the hopping integral (which sets the energy 
scale, so we take t ^ 1 throughout this paper), U is the 
on-site Coulomb repulsion, fi is the chemical potential 
controlling the fermion density, and i runs over the sites 
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of a d-dimensional lattice; for the time being, we consider 
hopping between nearest neighbours only, as denoted by 
(...). The operators Cj^. and respectively create and 
annihilate a fermion with spin a on the (single) orbital 
centred at i, while ni^- = Cj^C;^. The Hubbard model 
describes the competition between opposing tendencies of 
itinerancy (driven by the hopping term), and localization 
(driven by the on-site repulsion). For a half- filled band 
(one fermion per site), it can be shown 2] that in the limit 
of strong repulsion the Hubbard Hamiltonian reduces to 
that of an isotropic antiferromagnetic Heisenberg model 
with an exchange coupling J = 4£^ /U . 

If we let the coupling U assume negative values, one 
has the attractive Hubbard model. Physically, this lo- 
cal attraction can have its origin in the coupling of 
fermions to extended (through polaron formation) or lo- 
cal phonons (such as vibrational modes of chemical com- 
plexes) Q. Given that its weak coupling limit describes 
the usual BCS theory, and that the real-space pairing is 
more amenable to numerical calculations, this model has 
also been useful in elucidating several properties of both 
conventional and high-tcmperature superconductors 3J. 

The Hubbard model [even in its simplest form, Eq. 
is only exactly soluble in one dimension, through the 
Bethe ansatz; correlation functions, however, are not 
directly available. In higher dimensions one has to re- 
sort to approximation schemes, and numerical techniques 
such as Quantum Monte Carlo (QMC) simulations have 
proven to be crucial in extracting information about 
strongly correlated fermions. 

Since the first Monte Carlo method for classical sys- 
tems was devised in the early 1950's ^ H, 1^; several 
QMC algorithms have been proposed. Their classifica- 
tion is varied, and depends on which aspect one wishes 
to single out. For instance, they can be classified accord- 
ing to whether the degrees of freedom lie in the contin- 
uum or on a lattice; or whether it is a ground-state or 
a finite-temperature framework; or whether it is varia- 



tional or projective; or even according to some detailed 
aspect of their implementation, such as if an auxiliary 
field is introduced, or if a Green's function is constructed 
by the power method. Excellent broad overviews of these 
algorithms are available in the literature, such as Refs. 
P, Sl , so here we concentrate on the actual details of the 
grand-canonical formulation with auxiliary fields leading 
to fermionic determinants j^]- We will pay special atten- 
tion to the implementations and improvements achieved 
over the years 0, 0, 0, 0, |3| ■ We will also have in 
mind primarily the Hubbard model, but will not discuss 
at length the results obtained; instead, illustrative refer- 
ences will be given to guide the reader to more detailed 
analyses, and we apologize for the omitions of many rele- 
vant papers, which was dictated by the need to keep the 
discussion focused on this particular QMC implementa- 
tion, and not on the Hubbard (or any other) model. 

In line with the tutorial purpose of this review, we in- 
troduce in Sec. m the basic ingredients of Monte Carlo 
simulations, illustrated for 'classical' spins. In this way, 
we have the opportunity to draw the attention of the in- 
experienced reader to the importance of thorough data 
analyses, common to both classical and quantum sys- 
tems, before embarking into the more ellaborate quan- 
tum formalism. Preliminary manipulations, approxima- 
tions involved, and the natural appearance of Green's 
functions in the framework are then discussed in Sec. 
mil In Sec. IIVI we describe the updating process, as 
well as the wide range of average quantities available to 
probe different physical properties of the systems. We 
then address two of the main difficulties present in the 
simpler algorithm presented so far: the still unresolved 
'minus-sign problem' (Sec. 01, and the instabilities at 
low-temperatures, (Sec. IVIjl for which two successful so- 
lutions are discussed. Conclusions and some perspectives 
are then presented in Sec. IVIII 



II. MONTE CARLO SIMULATIONS FOR 
'CLASSICAL' SPINS 

Our aim is to calculate quantities such as the partition 
function (or the grand-partition function) and various 
averages, including correlation functions; these quantities 
are obtained by summing over all configurations of the 
system. For definiteness, let us consider the case of the 
Ising model, 



(iJ> 



(2) 



where J is the exchange coupling and tr? = ±1. 

For a lattice with Ns sites there are 2^" possible con- 
figurations in phase space. Clearly all these configura- 
tions are not equally important: recall that the prob- 
ability of occurrence of a given configuration, C, with 
energy E{C), is oc exp [— _E(C)/fcsT]. Therefore, from 
the numerical point of view, it is not efficient to waste 



time generating all configurations. Accordingly, the ba- 
sic Monte Carlo strategy consists of an importance sam- 
pling of configurations 5., 6]. This, in turn, can be im- 
plemented with the aid of the Metropolis algorithm 0], 
which generates, in succession, a chain of the most likely 
configurations (plus fluctuations; see below). Starting 
from a random spin configuration C = \ai, a2, ■ ■ ■ , tTNg), 
one can imagine a walker visiting every lattice site and 
attempting to flip the spin on that site. To see how 
this is done, let us suppose the walker is currently on 
site i, and call C" the configuration obtained from C 
by flipping aj. These configurations differ in energy by 
AE = E{C') - E{C) = 2Jai Y!-^ ctj, where the prime in 
the sum restricts j to nearest neighbour sites of i. The ra- 
tio between the corresponding Boltzmann factors is then 



P{C') 
PiC) 



cx.p{-AE/kBT). 



(3) 



Thus, if AE < 0, C" should be accepted as a new config- 
uration in the chain. On the other hand, if AE > the 
new configuration is less likely, but can still be accepted 
with probability r'; this possibility simulates the effect of 
fluctuations. Alternatively 0, one may adopt the heat- 
bath algorithm, in which the configuration C" is accepted 
with probability 



(4) 



1 



In both cases, the local character of the updating of the 
configuration should be noted; that is, the acceptance of 
the flip does not influence the state of all other spins. As 
we will see, this is not true for quantum systems. 

The walker then moves on, and attempts to flip the 
spin on the next site through one of the processes just 
described. After sweeping through the whole lattice, the 
walker goes back to the first site, and starts attempting 
to flip spins on all sites again. One should make sure that 
the walker sweeps through the lattice many times so that 
the system thermalizes at the given temperature; this 
warming-up stage takes typically between a few hundred 
and a few thousand sweeps, but it can be very slow in 
some cases. 

Once the system has warmed-up, we can start 'mea- 
suring' average values. Suppose that at the end of the 
C-th sweep we have stored a value for the quantity 
A. It would therefore seem natural that after Na such 
sweeps, an estimate for the thermodynamic average (A) 
should be given by 



(5) 



However, the ^<^'s are not independent random variables, 
since, by construction, the configurations in the chain 
are somewhat correlated. One can decrease correlations 



between measurements by forming a group of G averages, 
Ai, A2, ■ ■ ■ , Ag, leading to a final average 



(6) 



3=1 



Ideally, one should alternate Nn sweeps without taking 
measurements with Na sweeps in which measurements 
are taken. 

Once the Ag's can be considered as independent ran- 
dom variables, the central limit theorem 15] applies, and 
the statistical errors in the average of A are estimated as 



6 A = 



{A') - {AY^ 



G 



1/2 



(7) 



In addition, systematic errors must be taken into con- 
sideration, the most notable of which are finite-size ef- 
fects. Indeed, one should have in mind that usually there 
are two important length scales in these calculations: the 
linear size, L, and the correlation length of the otherwise 
infinite system, ^ ^ \T—Tc\~". A generic thermodynamic 
quantity, Xl{T), then scales asflfij 



(8) 



where a; is a critical exponent to be defined below, and 
f{z) is a scaling function with very specific behaviours in 
the limits z <C 1 and z ^ 1. In the former limit it should 
reflect the fact that the range of correlations is limited 
by the finite system size, and one must have 



Xl{T) ~ const. • L^/", for i < 



(9) 



so that an explicit L-dependence appears. On the other 
hand, when correlations do not detect the finiteness of 
the system, the scaling function must restore the usual 
size-independent form. 



Xl{T)^\T-T,\-\ forL»^; 



(10) 



this defines the critical exponent x. This finite-size scal- 
ing (FSS) theory can be used in the analysis of data, 
thus setting the extrapolation towards the thermody- 
namic limit on firmer grounds. 

As a final remark, we should mention for complete- 
ness that the actual implementation of the Monte Carlo 
method for classical spins can be optimized in several as- 
pects, such as using bit-strings to represent states [iTj . 
and broad histograms to collect data p^ . 



III. DETERMINANTAL QUANTUM MONTE 
CARLO: PRELIMINARIES 



The above discussion on the Ising model was tremen- 
dously simplified due to the fact that the eigenstates of 
the Hamiltonian are given as products over single particle 



states. Quantum effects manifest themselves in the fact 
that different terms in the Hamiltonian do not commute. 
In the case of the Hubbard model, the interaction and 
hopping terms do not commute with each other, and, in 
addition, hopping terms involving the same site also do 
not commute with each other. Consequently, the parti- 
cles lose their individuality since they are correlated in a 
fundamental way. 

One way to overcome these non-commuting aspects 
is to notice that the Hamiltonian contains terms which 
are bilinear in fermion operators (namely, the hopping 
and chemical potential terms) and a term with four op- 
erators (the interaction term). Terms in the former cat- 
egory can be trivially diagonalized, but not those in the 
latter. Therefore, when calculating the grand partition 
function. 



Z = Tr 



(11) 



where, as usual, Tr stands for a sum over all numbers of 
particles and over all site occupations, one must cast the 
quartic term into a bilinear form. 

To this end, we first separate the exponentials with 
the aid of the Suzuki- Trotter decomposition scheme^l^, 
which is based on the fact that 



^At{A+B) 



AtAAtB 



0[{^Tf[A.B]], (12) 



for A and B generic non-commuting operators. Calling K, 
and V, respectively the bilinear and quartic terms in the 
Hubbard Hamiltonian, wc introduce a small parameter 
Ar through (3 = MAr, and apply the Suzuki- Trotter 
formula as 



M 



O [(At)2c/] . (13) 



The analogy with the path integral formulation of 
Quantum Mechanics suggests that the above procedure 
amounts to the imaginary-time interval (0, (3) being dis- 
cretized into M slices separated by the interval Ar. 

In actual calculations at a fixed inverse temperature /3, 
we typically set Ar — yjQ.12b/U and choose M — /3/Ar. 
As evident from Eq. H13|) . the finiteness of Ar is also a 
source of systematic errors; these errors can be downsized 
by obtaining estimates for successively smaller values of 
Ar (thus increasing the number, M, of time slices) and 
extrapolating the results to Ar — + 0. Given that this 
complete procedure can be very time consuming, one 
should at least check that the results are not too sen- 
sitive to Ar by performing calculations for two values of 
Ar and comparing the outcomes. 

Having separated the exponentials, we can now work- 
out the quartic terms in V. We first recall a well known 
trick to change the exponential of the square of an oper- 
ator into an exponential of the operator itself, known as 
the Hubbard-Stratonovich (HS) transformation. 



'2t: 



dx e 



(14) 



at the price of introducing an auxihary degree of freedom 
(field) X, which couples linearly to the original operator 
A; this result is immediately obtained by completing the 
squares in the integrand. However, before a transforma- 
tion of this kind can be applied to the quartic term of the 
Hubbard Hamiltonian, squares of operators must appear 
in the argument of the exponentials. Since for fermions 
one has = = 0, 1 (here we omit lattice indices to 
simplify the notation), the following identities - in terms 
of the local magnetization, m = n-^ — n^, and of the local 
charge, n = + ni - will suit our purposes: 

1 2 1 / , N 

"T"! = ~ (-'^^*^) 

The following points are worth making about Eqs. (|15|l . 
in relation to the HS transformation: (i) an auxiliary field 
will be introduced for each squared operator appearing 
on the RHS's above; (ii) the auxiliary field will couple to 
the local magnetization and to the local charge when, 
respectively, Eqs. p5a|l and (|15b|l are used; (iii) Eqs. 
p5a|l and (|15bp are respectively used for the repulsive 
and for the attractive models. 

Instead of the continuous auxiliary field of Eq. H14|l . in 
simulations it is more convenient to work with discrete 
Ising variables, s = ±1 ilOi|. Inspired by Eq. (|14|l . using 
Eq. H15a|) . and taking into accoimt that = 1, it is 
straightforward to see that 

-UATn^n^ ^ ig-^" \^ -s\m ^ 

2 ^ 

s=±l 

= -J2Y[ e-(-^+^K, U>0, 

(16) 

I 



where the traces are over Ising fields and over fermion 
occupancies on every site, and the product from i = M 
to 1 simply reflects the fact that earlier 'times' appear to 
the right. The time-slice index £ appears through the HS 
field Si{£) in 

vr{i)^^X'Js^{e) + ^^^l-^y (20) 

which are the elements of the Ng x Ng diagonal matrix 
\l'^{t). One also needs the Ng x Ng hopping matrix K, 



where the grouping in the last equality factorizes the con- 
tributions from the two fermion spin channels, cr = +, — 
(respectively corresponding to a =t, i), and 



coshA = el^l^^/2. (17) 



In the attractive case, the coupling to the charge [Eq. 
(I15b|l ] is used in order to avoid a complex HS transfor- 
mation; we then get 



^\U\Arn,n, ^(,A+^)(rv-i)^ < Q, 

(18) 

with A still being given by Eq. (|17|l . 

The HS transformations then replace the on-site inter- 
action by a fluctuating fleld s coupled to the magnetiza- 
tion or to the charge, in the repulsive or attractive cases, 
respectively. As a consequence, the argument of the ex- 
ponential depends explicitly on a in the former case, but 
not in the latter. As we will see below, this very im- 
portant difference is responsible for the absence of the 
'minus-sign problem' in the attractive case. 

We now replace the on-site interaction on every site of 
the space-time lattice by Eqs. (fTB|) or (fTH|) . leading to the 
sought form in which only bilinear terms appear in the 
exponential. For the repulsive case we get 



(19) 

I 

with elements 



—t if i and j are nearest neighbours, 
otherwise 



For instance, in one-dimension, and with periodic bound- 



Z = 



ary conditions, one has an L x L matrix, 



/ -< 

-t -t 
-t 



-t \ 





(22) 



\-t ••• -t 0/ 

With bihnear forms in the exponential, the fermions 
can be traced out of Eq. I|19|l according to the devel- 
opment of Appendices El and El From Eq. I|B6(I , with 
the spin indices reintroduced, and with the identification 

g-Arh"(£) ^ g-ArKg-ArV"(£) ^^^^ 



Trndet[l + B5,BX,_i...B?] , (23) 

{s} ^ 



where we have defined 



(24) 



in which the dependence with the auxiliary Ising spins 
has not been explicitly written, but should be under- 
stood, since they come in through the matrix V'^(^). In- 
troducing 



0-(H)^l + BXfBXf_i.. 
we arrive, finally, at 



(25) 



Z = 



Tr detOT({s}) -detOHls}) 



(26) 



where the last equality defines an effective 'density ma- 
trix', p{{s}). 



We have then expressed the grand partition function 
as a sum over Ising spins of a product of determinants. 
If the quantity under the Tr were positive definite, it 
could be used as a Boltzmann weight to perform impor- 
tance sampling over the Ising configurations. However, 
the product of determinants can indeed be negative for 
some configurations, giving rise to the 'minus-sign prob- 
lem'; see Sec. |3 

In order to implement the above framework, the need 
for another approximation is evident from Eq. 124|l : one 
needs to evaluate the exponential of the hopping matrix, 
which, in the general case, is neither an analytically sim- 
ple operation, nor efficiently implemented numerically. 
By considering again a one-dimensional system, we see 
that different powers of K, as given by Eq. (|22|l . generate 
many different matrices. However, we can introduce a 
'checkerboard break-up' by writing 



K = K(,") + Ki''), 



(27) 



such that K^"^ involves hoppings between sites 1 and 2, 
3 and 4,. . ., while K^*") involves hoppings between sites 2 
and 3, 4 and 5,. . .. We now invoke the Suzuki- Trotter 
decomposition scheme to write 

which leads to systematic errors of the same order as 
before. With this choice of break-up, even powers of 
Ki"\ a = a,6, become multiples of the identity matrix, 
and we end up with a simple expression. 



-ArK : 



K^") sinh (Art) + 1 cosh (Art) , (29) 



which is very convenient for numerical calculations. The 
breakup of Eq. 1)281) can be generalized to three dimen- 
sions as follows 



^-AtK _ ^-ArK<"' -AtK<"' -ArK<"' -ArKl''' -ArK<''' -ArK<'' 



-0[{Arr 



(30) 



where the separation along each cartesian direction is 
similar to that for the one-dimensional case. 

And, finally, we must discuss the calculation of average 
values. For two operators A and B, their equal-'time' 
correlation function is 



{AB) ^ iTrT r 
^ {«} 



ABl[< 



-ArK -ArV"(£) 



(31) 



If we now define the fermion average - or Green's function 



for a given configuration of the HS fields as 



I 



-Tr 



ABl[. 



piM) 

the correlation function becomes 

{AB)^^Tl (AB)[,yp{{s}). 

{s} 



, (32) 



(33) 



At this point it should be stressed the important 
role played by the Green's functions in the simulations. 
Firstly, according to Eq. (|33)l . the average value of an 
operator is straighforwardly obtained by sampling the 



corresponding Green's function over the HS configura- 
tions weighted by p{{s}). Secondly, as it will become 
apparent in Sec. II VI the single particle Green's function, 
(Cj^cl plays a central role in the updating process 
itself. In Appendix IgI we obtain this quantity as the ele- 
ment ij of an Ns x matrix [3, 0| : 



per DO" 



(34) 



which is, again, in a form suitable for numerical cal- 
culations. And, thirdly, within the present approach 
the fermions only interact with the auxiliary fields, so 
that Wick's theorem[23 holds for a fixed HS configu- 
ration^ [Til IT^: the two-particle Green's functions are 
then readily given in terms of the single-particle ones as 



){s] 



+ i4,^^Ms}K4:){s}, (35) 



where the indices include spin, but since the "f and | spin 
channels are factorized [c.j. Eq. (|32|) ]. these fermion av- 
erages are zero if the spins are different. All average val- 
ues of interest are therefore calculated in terms of single- 
particle Green's functions. 

As will be seen, unequal-time correlation functions are 
also important. We define the operator a in the 'Heisen- 
berg picture' as 



a(e) = a{T) 



(36) 



IV. DETERMINANTAL QUANTUM MONTE 
CARLO: THE SIMULATIONS 

The simulations follow the lines of those for classical 
systems, except for both the unusual Boltzmann weight, 
and the fact that one sweeps through a space-time lat- 
tice. With the parameters of the Hamiltonian, U and 
/i, as well as the temperature, fixed from the outset, we 
begin by generating, say a random configuration, {s}, 
for the HS fields. Since the walker starts on the first 
time-slice, we use the definition, Eq. (jSU, to calculate 
the Green's function at £ = 1. As the walker sweeps the 
spatial lattice, it attempts to flip the HS spin at every 
one of the Ng points. 

At this point, it is convenient to picture the walker 
attempting to flip the HS (Ising) spin on a generic site, 
i of a generic time slice, I. If the spin is flipped, the 

t I 

matrices and change due to the element ii of the 
matrices V^(l') and \/'^{t), respectively, being affected; see 
Eqs. (|20|l and (|24|1 . The expression for the change in the 
matrix element, as Si(£) —s\{£), is 

dV.^il) = V^{1, -s) - V{^{£, s) = -2Aasi(£) <5ij, (41) 

which allows us to write the change in B^ as a matrix 
product, 

B,- -> [B,n' - B^A,-(i), (42) 
where the elements of the matrix A^(i) are 



so that the initial time is set to r = At with this dis- 
cretization, and a^{£) ^ [a(i?)]t. In Appendix IdI we show 
that the unequal-time Green's function, for £i > £2, is 
given by |0 



Gfj(4;4 



(Ci,(£l)c],(^2)){s} 



'(^2 + l)]y,(37) 



in which the Green's function matrix at the i-th time 
slice is defined as 



with 



.'^(^)^[1 + A'^(^)]- 



A°'(^) = Bf_j_Bf_2 ■ ■ ■ B^BIj 



(38) 



(39) 



The reader should notice the order in which the prod- 
ucts of B's are taken in Eqs. (jSH), JSZl, and (jSHJ; in Eq. 
H37|l . in particular, the product runs from £2 + 1 to ^1, 
and not cyclically as in Eq. H39|) . Also, for a given config- 
uration {s} of the HS spins, the equal-time Green's func- 
tions do display a time-slice dependence, as expressed by 
Eq. H38|l : they only become (approximately) equal after 
averaging over a large number of configurations. 

For future purposes, we also define 



9u 



[l-g]ij, and G?=^[l-G]ij. 



(40) 



(0 ifj^k, 
[AJ(i)]jk=il ifj = k^i, (43) 

|^p-2Aas>(£) ifj = k = i. 

Let us now call {s}' and {s}, the HS configurations 
in which all Ising spins are the same, except for those 
on site (i, which are opposite. We can then write the 
ratio of 'Boltzmann weights' as 



, _ PiM') _ detpT ({sjO-detO^ {{s}') 
' ~ " detOT({s}).detOi({s}) 



(44) 



where we have defined the ratio of fermion determinants 
as 



detO" {{s}') 
detO'^ ({s}) ■ 



(45) 



It is important to notice that one actually does not 
need to calculate determinants, since is given in terms 
of the Green's function: 



Ra 



_ det [l + A"(^)A^(i)] _ 
det[l-h A'^(/)] 

= det[(l + A-(W(i))g"W] = 

= det [l+(l-g'^(^))(A,-(i)-l)] = 

-l + (l-gS(^))(e-2^^^'W-l). (46) 



The last equality follows from the fact that 

r,-(i)^A,-(i)-i 



(47) 



is a matrix such that all elements are zero, except for the 
i-th position in the diagonal, which is 7/(1) = e"^'^'^'*''^^^ — 
1. With this simple form for Ra-, the heat-bath algorithm 
is then easily implemented, with the probability of accep- 
tance of the new configuration being given by Eq. Q), 
with Eqs. lO and (|^ . 

If the new configuration is accepted, the whole Green's 
function for the current time slice must be updated, not 
just its element ii; this is the non-local aspect of QMC 
simulations we referred to earlier. There are two ways of 
updating the Green's function. One can either compute 
the 'new' one from scratch, through Eq. (|38(l . or iterate 
the 'old' Green's function, by following along the lines 
that led to Eq. (^BJ, which yields 

r(^) = [i + (i-g^(^))r,-(i)]"'g'^(^). (48) 

An explicit form for g'^{£) is obtained, e.g., by first using 
an ansatz to calculate the inverse of the matrix in square 
brackets above, 

[1 + (1 - g^(£))r^(i)] '' = i + x{i- g^ii)) r^i) (49) 

where x is to be determined from the condition that the 
product of both matrices in Eq. H49|l is 1. Using the fact 
that r^(i) is very sparse, one arrives at 



1 

Rrr 



(50) 



with being given by Eq. 146|l . The element jk of the 
Green's function is then updated according to 



l + [l-.9uW]7/(i 



(51) 



Alternatively, one could arrive at the same result by solv- 
ing a Dyson's equation for g^i^{i) jB]. 

After the walker tries to flip the spin on the last site 
of the £-th. time slice, it moves on to the first site of 
the {£ + l)-th time slice. We therefore need the Green's 
function for the (£ -|- l)-th time slice, which, as before, 
can be calculated either from scratch, or iteratively from 
the Green's function for the £-th time slice. Indeed, by 
comparing [g'^{(+ with [g'^(£)]~"'^, as given by Eq. 

(|38|l ■ it is easy to see that 



g^{£ + l) = BJg^e)[B^Y 



(52) 



which can be used to compute the Green's function in 
the subsequent time slice. 

While the computation from scratch of the new Green's 
function takes ^ operations, the iterations [Eq. H51|l 
and (|5^ ] only take ~ operations. Hence, at first sight 
the latter form of updating should be used. However, 



rounding errors will be amplified after many iterations, 
thus deteriorating the Green's function calculated in this 
way; see Sec. I VII A solution of compromise is to iterate 
the Green's function while the walker sweeps all spatial 
sites of £ 10) time slices, and then to compute a new 
one from scratch when the {£+ l)-th time slice is reached. 
This 'refreshed' Green's function will then be used to 
start the iterations again. Clearly, one should check g for 
accuracy, by comparing the updated and refreshed ones 
at the {£+ l)-th time slice; if the accuracy is poor, £ must 
be decreased. 

For completeness, we now return to the situation de- 
scribed in the beginning of this Section, in which the 
walker attempts to flip the HS spins on every spatial site 
of the £ = 1 time slice. Whenever the flip is accepted, the 
Green's function is updated according to Eq. (|51|l . After 
sweeping all spatial sites, and before the walker moves on 
to the first site of the £ — 2 time slice, we use Eq. ()52|l to 
calculate g'^(2). The walker then sweeps the spatial sites 
of this time-slice, attempting to flip every spin; if the flip 
is accepted, the Green's function is updated according 
to Eq. 1)51(1 . As mentioned above, this procedure is re- 
peated for the next £ time slices. After sweeping the last 
one of these, a new Green's function is calculated from 
the deflnition, Eq. (|38|l . For the subsequent £ time slices 
the iterative computation of g is used, and so forth. 

Similarly to the classical case, one sweeps the space- 
time lattice several times (warm-up sweeps) before cal- 
culating average values. After warming up, one starts 
the 'measuring' stage, in which another advantage of 
the present approach manifests itself in the fact, already 
noted, that all average quantities are expressed in terms 
of the Green's functions. This brings us to the discussion 
of the main quantities used to probe the physical prop- 
erties of the system, and how they relate to the Green's 
functions. 

We start with the calculation of the occupation, n. It 
is given by 



MN. 



MN. 



(53) 



where the first equality illustrates that the ensemble av- 
erage {i.e., over both fermionic and HS fields) {n^{£) + 
nii{£)) is itself averaged over all MNg sites of the space- 
time lattice. In the second equality the fermion degrees 
of freedom have been integrated out (hence the Green's 
functions), and the average over HS fields, represented 
by double brackets, is to be performed by importance 
sampling over Na HS configurations [c.f. Eq. H33|l ]. We 
then have 



1 



(54) 



It should be reminded that since these are grand- 
canonical simulations, the chemical potential must be 
chosen a priori in order to yield the desired occupation. 
The chemical potential leading to half filling in the Hub- 
bard model (with nearest-neighbour hopping only) is ob- 
tained by imposing particle-hole symmetry 0, l2lj | : one 
gets n = 1 for jj. = U/2 for all T and Ng. Away from 
half-filling, one has to determine by trial and error, 
which must be done repeatedly, since the dependence of 
n with /i (for given Ng and U) itself varies with tempera- 
ture. Overall, the dependence of n with ^ indicates that 
the Hubbard model is an insulator at half filling in any 
dimension: indeed, the compressibility 



1 
V 



dP 



1 



(55) 



where V and P are the volume and pressure, respectively, 
vanishes at /i = f7/2, so that no particles can be added to 
the system; see, e.g., Ref. for explicit data for n(/i) 
in two dimensions. 

The magnetic properties are probed in several ways. 
Since the components of the magnetization operator are 



= ^TCii+c-iCiT' 



nil - nil, 



(56a) 
(56b) 
(56c) 



Wick's theorem [Eq. 1(23] allows one to obtain the equal- 
time spin-density correlation functions as 



and a similar expression for the yy correlations. It is 
understood that the g's are to be evaluated at the same 
time slice as the m's on the LHS. 

If spin rotational symmetry is preserved, as in the case 
of a singlet ground state, Eqs. (|57a|l and (|57b|l should 
yield the same result. This is indeed the case for average 
values, but, in practice, transverse (i.e., xx and yy) cor- 
relations are much less noisy than the longitudinal ones 
; this can be traced back to the discrete HS transfor- 
mation, which singles out the z component by coupling 
the auxiliary fields to the 's. It should also be pointed 
out that a possible ferromagnetic state (either saturated 
or just a non-singlet one) would break this symmetry. 

Setting i = j in Eqs. (|57|l leads to the local moment. 



{ml)^{{m,ff + {mlf + {mlf) 



(58) 



which measures the degree of itinerancy of the fermions. 
For instance, considering again the Hubbard model at 
half filling, {(rn'()^), v = x,y, z, decreases from 1, in the 
frozen-charge state {i.e., when U 00), to 1/2, in the 
metallic state {U = 0). 



One can collect the contributions to the spin-spin cor- 
relation functions from different sites, and calculate the 
equal-time magnetic structure factor as 



S{q) 



(Si • s. 



(59) 



where = /2, since h = 1. S{q) shows peaks at 
values of q corresponding to the dominant magnetic ar- 
rangements. For the Hubbard model in one dimension, 
S{q) has peaks at q = 2fci? = nn, corresponding to quasi- 
long range order {i.e., algebraic spatial decay of correla- 
tions) in the ground state; one refers to this as a state 
with a spin- density wave |24l| . In two dimensions, and at 
half filling, the peak located at q = (tt, tt) j2^ signals 
Neel order in the ground state; away from half-filling, 
the system is a paramagnet [ill l2a| with strong AFM 
correlations at incommensurate q's |22 |. 

The non-commutation aspects imply that the suscep- 
tibility is given by |2^ 



X(q) 



,«q-(i-j) 



dr (Si(T).Sj 



(60) 



In the discrete-time formulation, the integral becomes a 
sum over different time intervals, which is carried out 
with the aid of the unequal-time Green's functions. Fo- 
cusing on one of the components of the scalar product 
above, say z, we have 



mi{h)mi{£2)) - [(2-5-' " 1) ((3»(^i)55'(^2))) + 

a.(T' 



+S,,,{{G?^{h;i2)G?^{ei;e2)))\. (61) 

Similarly to the magnetic properties, one can also in- 
vestigate whether or not there is the formation of a charge 
density wave j2^l27l |. With rij = ni| +ni|, the equal-time 
correlation function is given by 



E 

a, a' 



[{{9^i9!;))+S..'{{9!ig?i))], (62) 



in terms of which the charge structure factor becomes 



niTlj). 



It is also useful to define a charge susceptibility as 



Af(q) 



/9 



dr (ni(T)nj), 



(63) 



(64) 



where, in the discrete time formulation, the above aver- 
age values are given in terms of the unequal-time Green's 
functions as 

{ni{ii)ni{h)) = [((5u(^i)55'(^2))) + 

+<5..'((GJi(£i;^2)Gfj(4;4)))l;(65) 



it should be noted that the uniform charge susceptibihty 
is proportional to the compressibility, defined by Eq. H55|l . 

For the one-dimensional Hubbard model, the presence 
of a charge-density wave (CDW) is signalled by a cusp 
in C(q) and by a divergence (as T — > 0) in Af{q). While 
the insulating character at half filling prevents a CDW 
from setting in, away from half filling there has been 
some debate on whether the position of the cusp is at 
q — 2fci?, as predicted by the Luttinger liquid theorvp^. 
or at (7 = Akp, as evidenced by QMC data and Lanczos 
diagonalizations j29l] . In two dimensions, and away from 
half filling, the situation is even less clear. 

We now discuss superconducting correlations. The 
imaginary-time singlet pair-field operator is a general- 
ization of the BCS gap function |3(TI |. 



Ac(r) 



1 



(66) 



where the subscript C, in the form factor, /f(k), labels 
the symmetry of the pairing state, following closely the 
notation for hydrogenic orbitals; for instance, in two di- 
mensions one has 

s-wave: /^(k) — 1, (67a) 

extended s-wave: /s*(k) — cosk^ +cosfcy, (67b) 

dx2-y2-wa.ve: fd^2_y2(^) ~ coskx — cosky, (67c) 

dxy-wave: fd^y{k) = sin fc^^ sin fcj, , (67d) 

and several other possibilities, including triplet pairing, 
as well as linear combinations of the above. Ignoring for 
the time being any r-dependence, and introducing the 
Fourier transform of the annihilation operators, Eq. I|66j) 
can be written as 



A, = 



(68) 



with the site-dependent pair-field operator being given 
by 



(69) 



where the sum runs over lattice sites neighbouring i, with 
range and relative phase determined by 



/C(a) 



E/C(k)e 



(70) 



For the d3;2_^2-wave, for instance. 



/d,2_„2 (a) = ((^a^x + <^a,-x " 5^,y - , (71) 

where x and y are unit vectors along the cartesian direc- 
tions. 

Similarly to the magnetization, simple averages of pair- 
field operators vanish identically in a finite-sized system. 



so one again turns to correlation functions. We define 
the equal-time pairing correlation function as 



PaiJ)-(Al(i)Aca)+H.c 



(72) 



whose spatial decay is sometimes analysed as a probe for 
the superconducting state [sJl. Its uniform {i.e., q = 0) 
Fourier transform, 

^C = ]^E(4«^c(j)+H-C.), (73) 



has also been used, together with finite-size scaling argu- 
ments, to probe superconducting pairing |2lLl32ll8^. As 
before, it is a straightfoward exercise to express the above 
averages in terms of HS averages of Green's functions. 

Restoring the r-dependence, another useful probe is 
the uniform and zero-frequency (w = 0) pairing suscepti- 
bility, 



-E 



f dr (At(i,r)AJj) 



H.c. 



(74) 



A considerable enhancement of 11,^ over its non- 
interacting counterpart may be taken as indicative of a 
superconducting state; see, e.g., Refs. |13 and H^- 
fairness, there has been no real consensus on which is 
the best numerical probe of a superconducting state. It 
may be argued, for instance, that should be compared 
with the corresponding nonvertex correlation function; 
see, e.g., Refs. p^ls^ . 

While the above probes for superconductivity assume 
a given symmetry for the pairing state, there is an alter- 
native procedure which is free from such assumption |37l | . 
The Kubo linear response to a vector potential yla;(q, w) 
is given by the x-component of the current density. 



0'2;(q,w)) 
where 



^e2[(-X,)-A,,(q,c^)]A,(q,c^), (75) 



(76) 



is the kinetic energy associated with the a:-oriented links, 
and 



A:rx(q, j^m) 



E 



/3 



dr (j,(i,T)j,(0,0))e''i-'e- 



(77) 

is the imaginary-time and space Fourier transform [ojm = 
2iTmp is the Matsubara frequency |23|) of the current 
density correlation function. 



A,,(i,T) = (j,(i,r)j,(0,0)). 



(78) 



with 



r\ c. 



(79) 



As an exercise, the reader should express the current 
density correlation function in terms of HS averages of 
Green's functions. 

The superconducting properties of the system are as- 
sociated with its long wavelength static response {i.e., 
q ^ 0,0; = 0), but with a careful distinction in the way 
in which the order of the limits qx (Longitudinal) 
and Qy —^ (Transverse) are taken: 



and 



A^ 



lim A2;j;((7a;, = 0, iuj„i = 0), 



(80) 



lim AxxiQx 



0,qy,iuj„i = 0). 



As a result of gauge invariance ^37], one should always 
have 



A^ = i-Kx) 



(82) 



which can be used to check the numerics. On the other 
hand, the superfluid weight, (which is proportional 
to the superfluid density, ps) turns out to be 37J 



D, 



= Ps = i-Kx) - A^, 



(83) 



so that a Meissner state is associated with A^ ^ . 
In practice, {—K^) is calculated independently from A-^, 
and one checks whether the latter quantity approaches 
the former as qy — > 0. Since the calculations are per- 
formed on finite-sized systems, one should examine the 
trend of the data as the number of sites is increased; see 
Ref. 113 ■ 

If the system is not a superconductor, the current- 
density correlation function still allows us to distinguish 
between a metallic and an insulating state. Indeed, the 
limit q = 0, a; ^ (note the order of the limits, which 
is opposite to that used in relation to the superfluid den- 
sity!) of the conductivity determines whether or not the 
ground state of the system has zero resistance js^l- The 
real part of the zero temperature conductivity can be 
written in the form (Jxx = D5{ijj)-\-<7''^^^{uj), where crj.'jf(ijj) 
is a regular function of w, and the Drude weight, I?, de- 
scribes the DC response. The latter can be estimated at 
low temperatures by 1^31 



D 



[{^Kx) - A,,(q - 0, iuj.„, ^ 0)] , (84) 



which amounts to extrapolating discrete non-zero uorn 
data towards — 0- This can lead to a subtle in- 
terplay between the energy scales involved; see Ref. 
for examples. 

The above-mentioned criteria to determine whether 
the system is metallic, insulating, or superconducting are 
summarized in Ta ble □ The reader is urged to consult 
References 1^3, IsE for illustrative discussions on the 
many aspects involved in the data analyses. 



Table I: Criteria to determine whether a system is a super- 
conductor, a metal, or an insulator, from the behaviour of the 
current-current correlation function. 



Nature of the state 


Ds [Eq. (IHSJ] 


D [Eq. lUl] 


Superconducting 






Metallic 







Insulating 









At this point one should already be convinced that a 
very wide range of probes is available to scrutinize most of 
the possible physical properties of a model system. How- 
ever, two technical problems must be dealt with, which 
will be analysed in turn: the minus-sign of the fermionic 
determinant, and the instabilities appearing at low tem- 
peratures. 



V. THE MINUS-SIGN PROBLEM 

As mentioned before, the product of fermionic deter- 
minants is not positive definite, apart from in a few cases. 
The best known instance is the repulsive Hubbard model 
at half filling. To see why this is so, we employ a particle- 
hole transformation. 



such that 



4 - {-^)\o 



Hi, 



-^{(7 "1(7 



it 



(85) 



(86) 



and consider the fermionic determinant at the symmetric 
point, /i = U/2. One then has 

detQT = Tr TT e"^"^^ ^ ^'^"'^ = 
{"} 

{"} V 

^ e-^S.^'^.WdetO^ (87) 
where the tildes denote hole variables. Therefore, 
det 0^ ■ det 0^ > 0, for n = 1, [/ > 0. 



For the attractive Hubbard model, the lack of a- 
dependence in the discrete HS transformation [see the 
comments below Eq. ^] leads to OT({s}) = OH{s}), 
so that the product of determinants is positive for all fill- 
ings. Similar arguments apply to show that the fermionic 
determinant is always positive for the Holstein model for 
electron-phonon interactions ji^, 0| . 

In other cases, the fermionic determinant becomes neg- 
ative for some configurations. In order to circumvent 
this problem, recall that the partition function can be 
written as a sum over configurations, c = {s}, of the 
'Bohzmann weight', p{c) = det OT({s}) det 0^({s}). If 
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Figure 1 : The average sign of the product of fermionic deter- 
minants as a function of band filling, for the Hubbard model 
with U = 4: (a) 4 X 4 square lattice, for inverse tempera- 
tures l3 = 6 (circles) . 8 (squares), and 10 (triangles); adapted 
from Refs. [gl and (b) 6 x 6 (circles) and 8x8 (squares) 
square lattice, for fixed inverse temperature, /3 = 6; adapted 
from Ref. ,13]. (c) 4 x 4 x 4 (circles) and 6x6x6 (squares) 
simple cubic lattice, for fixed inverse temperature, P — 7. 
Lines are guides to the eye in all cases. 

we write p(c) = s(c)|p(c)|, where s(c) — ±1 to keep track 
of the sign of p(c), the average of a quantitity A can be 
replaced by an average weighted by |p(c)| as follows 

^ Eeb(c)k(c)A(c)]/EJp(c)l _ 
Ecb(c)k(c)]/Ecb(c)l 

E,P'(c)[.(c)] - {s),, ' ^ ^ 

where p'(c) = |p(c)|. Therefore, if the absolute value of 
p(c) is used as the Boltzmann weight, one pays the price 
of having to divide the averages by the average sign of 
the product of fermionic determinants^ {sign) = (s)p'. 
Whenever this quantity is small, much longer runs (in 
comparison to the cases in which {sign) ~ 1) are neces- 
sary to compensate the strong fluctuations in {A) p. In- 
deed, from Eq. Q we can estimate that the runs need to 
be stretched by a factor on the order of {sign)^^ in order 
to obtain the same quality of data as for {sign) ~ 1. 

In Fig.^a), we show the behaviour of {sign) as a func- 
tion of band filling, for the Hubbard model on a 4 x 4 
square lattice with J7 = 4, and for three different tem- 
peratures. One sees that, away from n = 1, {sign) is only 
well behaved at certain fillings, corresponding to 'closed- 
shell' configurations; i.e., those such that the ground 
state of the otherwise free system is non-degenerate • 
For the case at hand, the special fillings are 2 and 10 
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Figure 2: The logarithm of the average sign of the product of 
fermionic determinants as a function of inverse temperature, 
for the Hubbard model on a 4 x 4 square lattice with n — 
0.625, and for different values of the Coulomb repulsion: U = 
4 (circles) and 8 (squares). Lines are fits through the data 
points. Adapted from Ref. [ioll . 

fermions on 4 x 4 sites, leading to n — 0.125 and 0.625, 
respectively. At any given non-special filling, {sign) dete- 
riorates steadily as the temperature decreases, rendering 
the simulations inpractical in some cases. Further, Fig. 
n^b) shows that for a given temperature, the dip in {sign) 
gets deeper as the system size increases. One should 
note, however, that the position of the global minimum 
is roughly independent of the system size, which, for fill- 
ings away from the dip, allows one to safely monitor size 
effects on the calculated properties. In this respect, the 
three-dimensional model is much less friendly, as shown 
in Fig. ^c): for 0.3 < n < 1, {sign) is never larger than 
0.5 at the same filling for both system sizes. 

It is also instructive to discuss the dependence of {sign) 
with temperature, keeping fixed both the system size and 
the band filling. In Fig. [3 we show \n{sign) vs. (3 for the 
4x4 lattice at the closed-shell filling n = 0.625. (Actu- 
ally, these data have been obtained in Ref. [4fl] by means 
of a ground state algorithm, but they follow a trend sim- 
ilar to those obtained from the determinantal algorithm 
discussed here.) As U increases, the sign deteriorates 
even for this special filling. For other fillings the aver- 
age sign also decreases with U, and confirms the general 
expectation j4Qj that 

{sign) - e-'^^^'^, (90) 

where 7 depends on n and U. While for a given n, the 
dependence of 7 on J7 is monotonic, for a given U, 7 is 
smaller at the special fillings than elsewhere. 

The fundamental question then is how to prevent, or 
at least to minimize, this minus-sign problem. While one 
could be tempted to attribute the presence of negative- 



weight configurations to the special choice of Hubbard- 
Stratonovich transformations (HST's) used, it has been 
argued that even the most general forms of HST's are 
unable to remove the sign problem. It therefore seems 
that the problem is of a more fundamental nature. 

In order to pinpoint the origin of the problem, let us 
change the notation slightly and write the partition func- 
tion as 



(91) 



where a generic HS configuration of the i-th time slice 
is now denoted by 5'^ = {si{£), S2(^), ■ • • sNs(^)}; the set 
S = {5'i, 5*2, ■• ■ Sm} defines a complete path in the space 
of HS fields. Instead of applying the HST to all time slices 
at once, one can apply it to each time slice in succession, 
and collect the contribution from a given HS configura- 
tion of the first £ time slices into the partial path A3\ , 



Vi{{Si,S2,--- ,Si}) 



Tr 



,(92) 



where there are M — i factors of i3 = e"^'^'^ which have 
not undergone a HST. Clearly, Vo = Z since it corre- 
sponds to the situation in which no HS fields have been 
introduced yet; it is therefore positive. When HS fields 
are introduced in the first time slice, a "shower" of 2^= 
different values of Vi emerges from Vq; see Fig. 13 Each 
subsequent introduction of HS fields leads to showers of 
2^= values of Ve emerging from each Vg-i- In Fig. 13 
we only follow two representative partial paths: the top 
one is always positive, while the bottom one hits the £- 
axis at some £o. In the latter case, subsequent showers 
lead to both negative and positive values of the partial 
paths. According to the framework discussed in Sec. IIVI 
the simulations are carried out after performing the HST 
on all sites of all time slices. In the present context, this 
amounts to sampling solely the intersection of all possible 
paths with the vertical line at £ = M; see Fig. O If one 
were able to sum over all HS configurations, we would 
find that the number of positive Pm's would exceed that 
of negative Vai^s by an amount which, at low temper- 
atures, would be exponentially small. Since in practice 
only a limited set of HS configurations is sampled, it is 
hardly surprising that we find instances in which con- 
figurations leading to negative weights outnumber those 
leading to positive weights. This perspective also helps 
us to understand why simply discarding those negative- 
weighted configurations is not correct: the overall contri- 
bution of the positive-weighted configurations would be 
overestimated in the ensemble averages. 

The analysis of partial paths is at the heart of a recent 
proposal [43 to solve the minus-sign problem. It is based 
on the fact that when a partial path touches the 1^ = 
axis, it leads to showers which, when summed over all 
subsequent HS fields, give a vanishing contribution; see 




Figure 3: Schematic behaviour of partial paths (see text) 
as a function of their length. Only two representative paths 
emerging from the "shower" at ^ = are followed: one (full 
curve) leads to a positive contribution when it reaches £ — M, 
whereas the other (dashed curve) reaches P = at some £o- 



Fig. El In other words, a replacement of all B's in Eq. 
dnai by Y.s,^+^■■■ Ua B^iSeo+i) ■ ■ ■ does not change the 
fact that Veg — 0. Therefore, if one is able to follow the 
'time' evolution of the partial paths, and discard those 
that vanish at < M, only positive-weight configura- 
tions end up contributing at £ ^ M . However, though 
very simple in principle, this programme is actually quite 
hard to implement due to the need to handle the B factors 
without using a HS transformation. Zhang |43j| suggested 
the use of trial B^s, and preliminary results seem encour- 
aging. Clearly, much more work is needed along these 
lines in order to fully assess the efficiency and robustness 
of the method. 

Likewise, other recent and interesting proposals to deal 
with the minus sign problem need to be thoroughly scru- 
tinized. In the meron-cluster approach 0|, HS fields 
are introduced in all sites as usual, but during the sam- 
pling process (1) the configurations are decomposed into 
clusters which can be flipped independently, and (2) a 
matching between positive- and negative-weighted con- 
figurations is sought; see Ref. for details, and Ref. E3| 
for another grouping strategy. Another approach [4y|, 
so far devised for the ground-state projector algorithm, 
consists of introducing a decision-making step to guide 
walkers to avoid configurations leading to zero weight; 
it would be worth examining whether the ideas behind 
this adaptive sampling approach could also be applied 
to the finite temperature algorithm. In this respect it 
should be noted that Zhang's approach is closely related 
to the Constrained Path Quantum Monte Carlo [i^. in 
which the ground-state energy and correlation functions 
are obtained by eliminating configurations giving rise to 
negative overlaps with a trial state. 

In summary, at the time of this writing, QMC sim- 
ulations are still plagued by the negative-sign problem. 



Many ideas to solve this problem have been tested over 
the years, and they require either some bias (through 
resorting to trial states) or quite intricate algorithms 
(which render simulations of moderately large systems 
at low temperatures prohibitive in terms of CPU time), 
or both. We hope these recent proposals can be imple- 
mented in an optimized way. 



VI. INSTABILITIES AT LOW TEMPERATURES 

When the framework discussed so far is implemented 
into actual simulations, one faces yet another problem, 
namely, the fact that the calculation of Green's functions 
becomes unstable at low temperatures. As mentioned in 
the paragraph below Eq. (|52|l , the Green's functions can 
be iterated for about £ ~ 10 time slices, after which they 
have to be calculated from scratch. However, as the tem- 
perature is lowered, i.e., for /3 > 4, ^ must be decreased 
due to large errors in the iterated Green's function (as 
compared with the one calculated from scratch). One 
soon reaches the situation in which the Green's func- 
tion has to be calculated from scratch at every time slice 
{i.e., £ = 1). It should be noted that this corruption 
only occurs as one iterates from one time slice to an- 
other, and not in the updating stage within a given time 
slice '13]. Worse still is the fact that as the temperature 
is lowered further, the Green's function cannot even be 
calculated from scratch, since = 1 + A"^ {£) becomes 
so ill-conditioned that it cannot be inverted by simple 
methods. For instance, in two dimensions and for U = 0, 
the eigenvalues of 0'^ range between ~ 1 and ~ e'^^; for 
U we therefore expect the ratio between the largest 
and smaller eigenvalues of 0'^ to grow exponentially with 
M , thus becoming singular at low temperatures. 

Having located the problem, two solutions have been 
proposed, which we discuss in turn. 



A. The Space-time Formulation 

The approach used so far can be termed as a space for- 
mulation, since the basic ingredient, the Green's function 
[or, equivalently, the matrices A*^ and 0*^, c.f. Eqs. (|26|l . 
()38|l . and (jSHJ] is an Ng x Ng matrix. However, within 
a field-theoretic framework, if space is discretized before 



integrating out the fermionic degrees of freedom 0, the 
matrices 0*^ are blown up to 
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(93) 



which is an (NgM) x (NgM) matrix (since each of the 
M X M entries is itself an Ng x Ng matrix) ; one still has 



i 



(94) 



det 0'^ = det [1 + BIj • • • BJ] , 
as in Eq. (|23J), and 

Z= (^-j Trdet6T(M).det6HW)- (95) 

Taking the inverse of 0°' yields immediately the space- 
time Green's function matrix, 



0° 



(96) 



This space-time formulation has the advantage of 
shrinking the range of eigenvalues of 0*^: the ratio be- 
tween its largest and smallest eigenvalues now grows 
only linearly with M, thus becoming numerically sta- 
ble Though this approach has been extremely use- 
ful in the study of magnetic impurities i it ^o^^ slow 
down the simulations when applied to the usual Hubbard 
model "31. Indeed, dealing with {NgM)x{NgM) Green's 
functions would require {NgM)"^ operations per update, 
N^M'^ operations per time slice, and, finally, {NgM)^ 
operations per sweep of the space-time lattice. Sweeping 
through the space-time lattice with g'^ is then a factor of 
slower than sweeping with g'^. 

A solution of compromise between these two formula- 
tions was proposed by Hirsch ■ Instead of using one 
time slice as a new entry in 0°' [Eq. (|93|l ]. we collapse 
Mo < M time slices into a new entry. That is, by tak- 
ing Mq = M/p, with p an integer, one now deals with 
{Ngp) X (Ngp) matrices of the form 



0^(1) 







B 



2Mo 



B 



Mo + l 



V 

in terms of which the partition function is calculated as 



^pMo 



B 



(p-l)Mo + 



A 



(97) 



in Eq. (|95|l . The label 1 of 0^^^ indicates that the product 



of B's start at the first time slice of each of the p groupings. As a consequence, the time-dependent Green's function 
sub-matrices, G"^ {£i; £2), connecting the first time sHce of each grouping with either itself or with the first time shce 
of subsequent groupings, are readily given by 121 



gX/o(i) 



/ G(l,l) 
G(A/o + l,l) 



G(l,Mo + l) 
G(Mo-l-l,Mo-f 1) 



G(l,(p-l)Afo + l) \ 
G{Mo + l,{p-l)Mo + l) 



VG((p-l)Afo + l,l) G((p-l)Mo + l,Afo + l) ••• G((p-l)Afo + l,(p-l)Afo + l)/ 

I 



where the spin indices in G have been omitted for simplic- 
ity. The Green's functions connecting the ^-th time slices 
of the p groupings are similarly obtained from the inver- 
sion of Omo I which, in turn, is obtained by increasing 
all time indices of the B's in Eq. l|i?7| by £ — 1. 

In the course of simulations, one starts with the 
time-expanded Green's function calculated from scratch 
through Eq. H98|l . and sweeps all sites in the first time 
slice of each of the p groupings. Each time the move is 
accepted, the Green's function at that time slice is up- 
dated according to Eq. if^ . After sweeping over these 
sets of lattice sites, one iterates the Green's functions to 
obtain the elements of gMo (2) through 0| 

Q{ti + 1,^2 + 1) = B,,G(4,4)B7/. (99) 

We then sweep through all sites of the second time slice 
of each of the p groupings, and so forth, until all time 
slices are swept. Note that since Eq. is only used 
A/q times (as opposed to M times) it does not lead to 
instabilities for Afo small enough. 

Each Green's function update with Eq. (|51|l requires 
~ {NspY operations; sweeping through all spatial sites 
of one time slice in each of every p groupings therefore 
requires ~ {Nsp)^ operations. Since there are Mq slices 
on each grouping, one has, finally, a total of ~ N^Mp^ 
operations per sweep, which sets the scale of computer 
time; this should be compared with N^M for the original 
implementation, and {NgM)'^ for the impurity algorithm. 
The strategy is then to keep A/q ^ 20 and let p increase as 
the temperature is lowered. With this algorithm, values 
of /3 ~ 20 — 30 have been achieved in many studies of 
the Hubbard model; see, e.g., Refs. [H IH El S HI 
113, 113, 1^ • It should also be mentioned that since 
the unequal-time Green's function is calculated at each 
step, this space-time formulation is especially convenient 
when one needs frequency-dependent quantities |37| . 



B. Matrix-decomposition stabilization 

Let us assume that Mq of the B matrices can be mul- 
tiplied without deteriorating the accuracy. One then de- 
fines il^i3 



which, by Gram-Schmidt orthogonalization, can be de- 
composed into 

A?(£) = U^D^^R^, (101) 

where UJ' is a well-conditioned orthogonal matrix, is 
a diagonal matrix with a wide spectrum (in the sense dis- 
cussed at the beginning of the Section), and is a right 
triangular matrix which turns out to be well-conditioned 

Using the fact that is well-conditioned, we can mul- 
tiply it by another group of ATq matrices, 

Q = ^l+2Mo^1+2Mo-l ■ ■ ■ Bf+Mo + lUl, (102) 

without compromising the accuracy. We then form the 
product 

Q' = QD^, (103) 

which amounts to a simple rescaling of the columns, with- 
out jeopardizing the stability for the subsequent decom- 
position as 

Q' = U^D^R^, (104) 

where the matrices U2 , D2 , and Rj satisfy the same con- 
ditions as in the first step, Eq. ifTUTIl . With R| = R^Rf , 
we conclude the second step with 

A-(£) = U^D^R-. (105) 

This process is then repeated p = M/Mq times, to finally 
recover Eq. 139|l . in the form 13] 

A-(^) = a; = U^D-R;. (106) 

The equal-time Green's function is therefore calculated 
through Eq. as before, but care must be taken when 
adding the identity matrix to A, since the latter involves 
the wide-spectrum matrix Dp. One therefore writes 

1 = U^[U^]"'r^[R^]-\ (107) 
so that the inverse of Eq. H38|l becomes 



A^(£) _ B'^^^j^B'^_^_^,f^^_^ • ■ ■ , 



(100) 



[g'^(^)]-^ = U^P-R^, 



(108) 



with 



Acknowledgments 



(109) 



We now apply another decomposition to P, the result of 
which is multiplied to the left by Up, and to the right by 
R^, in order to express the Green's function in the form 




(110) 



In the course of simulations, one proceeds exactly as 
in Sec. IIVI by updating the Green's function through 
iterations, which takes up ~ N^M operations. Again, 
the iteration from one time slice to another is limited to 
about I time slices, and it turns out that a significant 
fraction of the computer time is spent in calculating the 
Green's function from scratch. Indeed, a fresh g is cal- 
culated M / H. times, each of which involving p decomposi- 
tions, each of which taking about N'^ operations. There- 
fore, taking i ^ p, the dominant scale of computer time 
is ~ N^p^ instead, which is about a factor of M faster 
than the space-time algorithm, though without the bonus 
of the unequal-time Green's functions. This matrix de- 
composition algorithm has also been efficiently applied 
to several studies of the Hubbard model, and values of 
^ 20 — 30 have been easil y ac hieved; see, e.g., Refs. 



p ^ ZU — M nave been easil y ac 
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Appendix A: GROUPING PRODUCTS OF 
EXPONENTIALS 

One frequently needs to group products of exponen- 
tials into a single exponential. Here we establish a crucial 
result for our current purposes: the product of bilinear 
forms can be grouped into a single exponential of a bi- 
linear form. 

To show this result , let the generic Hamiltonian be 
expressed as (we omit spin indices) 



(Al) 



where h is a matrix representation of the operator Ti. 
The 'time' evolution of Heisenbcrg operators [c.f., Eq. 
(jSnH , satisfies a differential equation, whose solution can 
be found to be 



4(r) = E[^""].,4. 



(A2) 



VII. CONCLUSIONS 

We have reviewed the Determinantal Quantum Monte 
Carlo technique for fermionic problems. Since the semi- 
nal proposal by Blankenbecler, Scalapino, and Sugar 
over twenty years ago, this method has evolved tremen- 
dously. Stabilization techniques allowed the calculation 
of a variety of quantities at very low temperatures, but 
the minus-sign problem still plagues the simulations, re- 
stricting a complete analysis over a wide range of band 
fillings and coupling constants. In this respect, it should 
be mentioned that other implementations of QMC, such 
as ground-state algorithms (see Ref. 8]), also suffer from 
this minus-sign problem. An efficient solution to this 
problem would be a major contribution to the field. 

Most of the discussion was centred in the simple Hub- 
bard model, but advances have been made in the QMC 
study of other models, such as the extended Hubbard 
model the Anderson impurity model ji^l, and the 
Kondo lattice model The first applications of QMC 
simulations to disordered systems have appeared only re- 
cently; see, e.g., Refs. '38*, 3?|. With the ever increasing 
power of personal computers and workstations, one can 
foresee that many properties of these and of more ellab- 
orate and realistic models will soon be elucidated. 



If we now take the linear combination 

4 = Ei3„cl, (A3) 

i 

we see that 

i j 
i 

where the last passage has been used as the definition of 
B. 

Now we introduce the product over time slices, through 



(A5) 



so that the h's now acquire i'-labels, and we have 



- EE 



(A6) 



2 ] 



where, again, the last passage has been used to define h. 
We can then write 



Appendix C: EQUAL-TIME GREEN'S 
FUNCTIONS 



(A7) 



so that U indeed behaves as a single exponential of the 
effective bilinear Hamiltonian 7i, whose elements are de- 
fined by 



-ATh(£) 



-Arh 



Appendix B: TRACING OUT EXPONENTIALS 
OF FERMION OPERATORS 

We now use the results of App. El to show that the 
trace over the fermionic degrees of freedom can be ex- 
pressed in the form of a determinant. Given that 7i is 
bilinear in fermio n opera tors, it can be straightforwardly 
diagonalized 



rmio n op era 

HiEflil 



ii- — ^ ^ ^aC^C^, 



(Bl) 



where, for completeness, the relations between 'old' and 
'new' fermion operators are 



a\i)Ci, 



and 



i 

the inverse relations are 

Ct = y^(^|a)ca, 

a 

and 

4 = 



(B2) 



(B3) 



(B4) 



(B5) 



Given the diagonal form, the fermion trace is then im- 
mediately given by 

e 

a 

= 11(1 + ^-^^^°) = 

a 

= det [l + e-^^^' 



= det 



1+n^ 



(B6) 



The equal-time Green's function, defined as the 
fermion average for a given configuration of the HS fields, 
is (we omit the fermion spin label) 



1 



(CI) 



(A8) where one assumes, for the time being, that and 
carry no time label. 

As in the calculation of the partition function, Eq. 
(|19|l . we write exp(— /37i) as a product of M factors 
exp(— ArTi). We then introduce the matrix represen- 
tation, 



to write 



TrDM---Dic^c] 



TrD 



M 



(C2) 



(C3) 



In order to evaluate Tr , we first change to the diagonal 
representation through Eqs. (jB4p and (|B5p . to write 



TrD 



M ■ 



When we take the Tr on the diagonal basis, c^/cj^,,, act 
on the bra to the left, so a nonzero contribution only 
occurs for a' = a". Then, all states but jo;') contribute to 
the product with (1 -I- exp(— Area,)) [see Eq. IjBGIl ] which 
cancels with the same term in the denominator. We are 
then left with 

= [gh- (C5) 



1 



1 + Bm • • • Bl 



We can generalize the above result and write the 
Green's function at the ^-th time slice as 



g"(£)= [l + B,_iB,_2 



BiB 



M 



(C6) 



The reader should convince himself/herself that the av- 
erage of g(^) over all HS configurations should be inde- 
pendent of £. 



Appendix D: UNEQUAL-TIME GREEN'S 
FUNCTIONS 



The unequal-time Green's functions are defined as (we 
omit the spin indices, for simplicity) 



G,,{h;£,) ^ (c,(4)c](4)){4 = lrre-'^^e^^^^^c,e-^^^^«e^^^^«cte-^^^-«. (Dl) 

As in the calculation ol the partition lunction, we write the first exponential to the right of the Tr as a product of 
M factors e^'^'^^, so that there are, in effect, M — di of these factors before q; by the same token, there are £i — £2 
such factors between c^ and , and ^2 factors to the right of c| . We then use the cyclic property of the Tr to bring 
these latter factors to the front, and introduce — I2 factors of e~^'^'^ followed by the same number of e^'^'^ factors 
before The Suzuki- Trotter decomposition is then used to break the kinetic energy and interaction terms, and the 
discrete HS transformation is applied to the exponentials involving the interaction; this adds a time-slice label to the 
resulting group of exponentials, and we can write 

G^J{ll■,^2) = ^TrDe,---DiDM---De+iDf-De,+i [D,, • • • [D,, • • • c], (D2) 

where -D^ is given by Eq. The product of D's can be expressed in its diagonal representation, 

=eS.I«>-"("l, (D3) 
which can also be used to express Ci, as in Eq. (|B4|I . Since 



e '^''C'Y if a = 7, 



* if a 7^ 7, 



(D4) 



we have 

J2 m [De. ■ ■ ■ De,+iV' [D,, ■ ■ ■ Dt,+,] = ^ (*|7) e-^-^^-^-^c^. (D5) 

7 7 

Thus, since 1 — ?1q = 0, 1, we can write 

^-e,(l-n,) ^ (^|B,^...B,, + i|7), (D6) 

which leads to 

[D,, ■ ■ ■ Di,+,r' c, [Di, ■ ■ ■ Di,+i] = ■ ■ ■ B^,+i|fc) Cfc, (D7) 

k 

where the sum is in site indices. 
With linil), Eq. ijnil) becomes 

(^i;^2) - ^ ^ (i|Bf, • • • Be,+i\k) TrDe, ■ ■ ■ D^Dm ■ ■■Dt,+ic^c], (D8) 
fc 



which can, finally, be expressed in terms of the equal-time 
Green's functions, Eq. IjCfip . as 



; + l 



;(^2 + 1)], 



(D9) 



The reader should check that an analogous calculation, 
still for li > ^2, leads to 

G.,(^i;^2) ^ (cl(fi)c^.(^2)>M - 

[l-g(^2 + l)] (B,,B,,_i...B,,+i)-i . 

(DIO) 



In the main text, the Green's functions given in Eqs. 
D9|l and (|D10|I acquire a spin index. 
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